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I  INTRODUCTION 


A.  BACKGROUND 

The  detection  of  a  known  deterministic  signal  in  additive  Gaussian  noise  is 
solved  using  a  correlator  or  matched  filter  [45,  55,  62,  12  Part  I].  For  the  case 
of  a  random  signal  in  white  noise,  a  two  stage  process  is  utilized  to  first  estimate 
the  signal  and  then  correlate  this  estimate  with  the  received  waveform.  This 
detection  scheme  is  referred  to  as  the  estimator-correlator  [47-50,  53,  12  Part 
III].  The  advantage  of  this  later  approach  is  its  practical  implementation 
structure;  ie,  the  best  obtainable  estimate  of  the  signal  is  used  in  the  correlation 
stage.  As  noted  by  Kailath  [54],  however,  the  method  described  in  [50]  fails  when 
the  signal  and  noise  are  correlated  and  when  signals  have  nonzero  means.  These 
limitations  were  overcome  for  continuous-time  processes  by  Kailath  [53,  56] 
using  the  innovations  representation.  In  addition,  the  approach  applied  to  non- 
Gaussian  desired  signal  processes.  This  detection  problem  is  expressed  as 

Hi:  x(t)  =  s(t)  +  w(t)  (1-la) 

Ho:  x(t)  *  w(t)  (1-lb) 

where  Hi  i=0,l  are  the  hypotheses  for  signal  absent  and  present,  respectively;  s(t) 
is  the  (not  necessarily  Gaussian)  random  signal  and  w(t)  is  Gaussian  white  noise. 
Under  Hi,  the  optimal  MMSE  causal  estimate  of  s(t)  is 

s(t)  »  E  [s(t)|  x(t"),  t'  <  t]  (1-2) 

Using  the  definition  of  the  innovations  [57];  ie. 

v(t)  =  x(t)  -  s(t),  (1-3) 

the  innovations  theorem  for  continuous-time  processes  (ie,  the  covariance  of  v(t) 
and  w(t)  are  identical),  and  the  property  that  v(t)  is  Gaussian  if  w(t)  is  Gaussian 
[53],  Kailath  transformed  the  detection  problem  of  eq  (1-1)  to  the  equivalent 
detection  problem: 
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Hi:  x(t)  -  £(t)  +  v(t) 
Ho:  x(t)  *  v(t) 


(Ma) 


(Mb) 

In  this  form  of  the  detection  problem,  s(t)  is  viewed  as  a  conditionally-known 
signal  and  the  likelihood  ratio  becomes  [53] 

LR  =  exp  [  £  £(t)  x(t)  dt  - 1/2  J  tHfy  dt].  (1-5) 

This  likelihood  ratio  has  the  same  form  as  the  estimator-correlator  for  the 
Gaussian  signal  in  the  white,  Gaussian  noise  case  except  that  the  first  integral  in 
eq(l-5)  is  an  Ito  integral.  Eq(l-5)  is  significant  because  it  shows  the  structure  of 
the  non-Gaussian  detection  problem.  In  addition,  although  eq(l-5)  was  derived  in 
[53]  assuming  s(t)  and  w(t)  as  statistically  independent,  this  restriction  was  later 
[56,58]  relaxed  to  "future  w(*)  independent  of  past  s(0";  ie.  the  signal  could  be 
correlated  with  past  noise.  It  is  interesting  to  note  that  die  innovations  approach 
does  not  use  a  Karhunen-Loeve  expansion  which  requires  s(t)  and  w(t)  to  be 
statistically  independent.  This  later  generalization  indicates  that  the  innovations 
approach  may  offer  a  robustness  in  detection  problems  involving  feedback, 
multipath,  etc. 

Unlike  the  continuous-time  case,  discrete  innovations  processes  do  not 
retain  the  property  of  Gaussianity  unless  both  w(t)  and  s(t)  are  Gaussian.  In 
addition,  the  covariance  of  the  innovations  is  not  equivalent  to  that  of  the  white 
noise  [53,57].  Thus,  the  likelihood  ratio  for  the  discrete  case  does  not  have  an 
estimator  correlator  structure  to  provide  a  test  statistic.  For  this  case,  a 
likelihood  ratio  using  discrete-time  innovations  processes  was  developed  [4,  6]. 
In  [4],  it  is  shown  that  for  non-Gaussian  processes,  this  likelihood  ratio,  called  the 
Innovations-Based  Detection  Algorithm  (IBDA),  is  a  close  approximation  to  the 
likelihood  ratio  for  this  detection  problem.  Furthermore,  the  structure  of  the 
likelihood  ratio  is  obtained  by  assuming  a  specific  parametric  model  of  the 
observation  processes  for  each  hypothesis  and  designing  a  prediction  error 
filtering  stage  for  each.  Adaptive  algorithms  are  then  used  to  estimate  the 
unknown  model  parameters.  Since  the  order  of  the  model  under  Hi  is  greater 
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than  that  under  Ho  (and  with  possibly  distinct  coefficients)  the  error  output  from 
the  filter  designed  for  Ho  given  that  Hi  is  true  will  be  greater  than  the  error 
output  from  the  filter  designed  for  Hi.  These  error  differences  are  used  in  the 
likelihood  ratio  to  raise  or  lower  its  value  relative  to  the  threshold  level.  In  [59], 
it  is  shown  that  this  likelihood  ratio  is  a  generalization  of  several  important 
special  cases.  These  include  (1)  the  detection  of  a  deterministic  signal  in  additive 
white  Gaussian  noise,  (2)  the  detection  of  a  non-white  Gaussian  signal  process  in 
additive  white  Gaussian  noise,  (3)  the  moving  target  detection  (MTD)  algorithm 
for  the  detection  of  a  deterministic  signal  in  non- white  Gaussian  noise  [60,  61] 
and  (4)  the  detection  of  a  non-white  Gaussian  signal  process  in  additive  non-white 
Gaussian  noise  [62],  It  is  also  shown  [59]  that  the  IBDA  contains  the  algorithm 
developed  [63]  for  the  detection  of  a  deterministic  signal  in  additive  non-white 
Gaussian  noise  of  unknown  correlation  statistics.  An  experimental  investigation 
using  two  implementations  of  the  IBDA  was  conducted  in  [8]  for  an  airport 
surveillance  radar  application  and  performance  comparisons  made  with  three 
MTD  algorithms. 

Several  early  analyses  involving  model-based  parametric  detection 
approaches  are  presented  in  [40,  44,  62,  64]  with  special  notice  given  to  [51]. 
Most  of  die  interest,  however,  appears  to  have  taken  place  within  the  past  few 
years  [4,  6,  8,  14,  21,  23,  41,  43,  45,  63].  The  principal  advantage  of 
characterizing  the  observation  processes  for  each  hypothesis  via  a  parametric 
model  is  that  well  known  algorithms  can  be  utilized  to  estimate  the  parameters. 
In  [45],  a  likelihood  ratio  test  was  considered  for  two  known  autoregressive  (AR) 
models.  In  [41],  a  more  general  formulation  considers  AR  and  autoregressive- 
moving  average  (ARMA)  models  to  detect  a  Gaussian  signal  in  white  Gaussian 
noise,  both  with  unknown  statistics.  Thus,  the  IBDA  [4,  6]  mentioned  above 
considers  the  more  general  problem  of  detection  on  a  non-Gaussian  signal  in  non¬ 
white  plus  white  noise.  As  noted  in  [16],  however,  different  modeling  approaches 
will  generally  yield  differences  in  receiver  performance.  In  addition,  the 
problem  of  modeling  observation  processes  in  this  detection  scheme  is  more 
difficult  than  model-fitting  a  time  series  process.  This  is  due  to  the  fact  that  in 
the  detection  problem,  one  set  of  observation  data  is  given  and  the  problem  is  to 
determine  which  of  the  two  filters  is  estimating  the  parameters  properly.  In  [14, 
16],  Zhang  investigated  the  detection  performance  improvement  over  the  IBDA 
[6]  in  a  radar  application,  when  ’a  priori'  information  was  used  to  predetermine 


the  process  parameters;  specifically,  a  reference  channel  which  provided  data 
from  range  cells  adjacent  to  the  "test  cell"  was  used  to  determine  the  filter 
coefficients  under  Ho.  A  significant  detection  performance  improvement  over 
the  IBDA  was  reported.  It  thus  appears  that  performance  improvements  can  be 
made  in  the  model-based  detection  schemes  through  investigations  of  alternative 
algorithms  and  implementation  schemes  for  the  detection  problem  with  unknown 
statistics. 

In  the  development  of  the  model-based  detection  approaches,  the  emphasis 
is  placed  on  characterizing  the  observation  data  received  under  each  hypothesis 
with  approximate  models.  If  the  models  fail  to  fit  the  physical  processes, 
performance  degradations  will  result.  Model  fitting  of  observation  data  via  time 
series  analysis  has  received  considerable  attention  [19,  65,  66,  67,  68,  69,  71,  72]. 
The  emphasis  in  the  analysis  proposed  here  will  be  the  development  of 
multichannel  model-based  detection  approaches.  Multichannel  time  series 
analyses  have  also  been  investigated  [28,  35,  66  Part  II,  67,  69,  71,  72]  with 
emphasis  in  the  areas  of  geophysics,  biophysics  and  economics.  Furthermore, 
multichannel  algorithms  for  parameter  estimation  have  received  some  attention 
[3,  9,  24,  25,  31,  70,  73].  The  prime  consideration  for  this  analysis  will  be  the 
investigation  of  the  potential  for  improved  detection  performance  of  multichannel 
model-based  detection.  Thus,  the  applicability  of  multichannel  time  series  models 
and  the  performance  of  multichannel  parameter  estimation  algorithms  to  the 
detection  problem  proposed  here  is  essential. 
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B.  RESEARCH  OBJECTIVE 


The  proposed  investigation  will  consider  the  binary  multichannel  detection 
problem  for  an  unknown  random  signal  vector  in  additive  non-white  interference 
plus  white  Gaussian  noise.  The  observation  processes  will  be  assumed  to  have  an 
arbitrary  correlation  in  time  and  space  (ie,  across  the  channels).  Initially, 
Gaussian  random  processes  will  be  considered,  but  a  generalization  to  non- 
Gaussian  processes  will  be  developed. 

The  principal  research  objective  is  the  investigation  of  multichannel  model 
based  detection  methods  utilizing  estimation  to  determine  their  performance 
relative  to  the  single  channel  case.  In  this  approach  to  the  detection  problem,  it  is 
assumed  that  the  underlying  physical  mechanisms  which  give  rise  to  the  observed 
random  processes  obtained  under  each  hypothesis  can  be  represented  by  a 
mathematical  model  which  approximates  its  statistical  characteristics.  We 
therefore  make  a  distinction  between  the  model  of  the  process  (i.e.,  synthesis)  and 
the  estimation  process  (analysis).  Specifically,  we  will  give  prime  consideration 
to  multichannel  autoregressive  processes  [9,  24,  25,  28-31]  as  the  process  model 
description  in  the  proposed  investigation.  The  consideration  of  the  multichannel 
approach  is  based  on  the  contention  that  the  coefficients  of  the  AR  processes  are 
distinct  for  each  of  the  two  hypotheses  (i.e.,  signal  present  or  absent).  The 
approach  used  in  the  model  based  detection  method  is  the  selection  between  the 
hypothesis  based  upon  measures  which  are  sensitive  to  the  differences  in  the 
process  coefficients  for  each  hypothesis.  Therefore,  a  likelihood  ratio  is  sought 
which  is  sensitive  to  differences  between  the  parameters  of  the  processes.  In  this 
context,  we  view  the  multichannel  processes  arising  from  physical  mechanisms 
such  as  those  which  may  yield  additional  (although  partially  correlated) 
information  about  the  processes.  From^the  model  based  approach,  this  new 
information  can  be  utilized  to  provide  a  better  distinction  between  the  process 
parameters  under  each  hypothesis.  The  extraction  of  this  new  information  is 
achieved  in  the  processing  of  the  observation  data  to  remove  the  redundant  (i.e. 
correlated)  information.  This  is  achieved  via  estimation  methods  which  "whiten" 
the  data  in  time  and  space  (i.e.,  channels).  These  uncorrelated  processes  contain 
all  the  useful  information  about  the  processes  in  a  compact  form  and  are  utilized 
to  determine  a  sufficient  statistic  for  the  hypothesis  determination;  i.e.,  a 
likelihood  ratio  can  be  developed  in  terms  of  these  transformed  processes  since 
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they  are  obtained  via  a  causal  and  causally  invertible  transformation  of  the 
original  observation  data.  As  shown  in  sections  IV  and  V,  these  processes  enable 
one  to  efficiently  perform  a  calculation  of  the  likelihood  ratio.  Furthermore,  the 
likelihood  ratio  considered  here  has  an  implementation  framework  amenable  to 
adaptive  processing  methods. 

Examples  of  application  areas  can  be  found  in  radar  technology, 
biomedical  applications,  geophysical  research,  image  processing,  data 
compression,  speech  analysis  and  channel  equalization.  In  the  radar  application, 
we  may  consider  the  processing  of  multi-sensor  pre-detection  data  or  dual- 
polarization  data  as  the  multichannel  processes.  A  significant  contribution  of  the 
detection  algorithm  considered  here  is  the  capability  to  utilize  signal  processing 
procedures  to  deal  with  partially  correlated  observation  data.  For  the  active 
radar  case,  we  view  the  multichannel  observation  data  as  processes  which  arise 
through  a  simultaneous  excitation  of  the  surveillance  volume  with  a  multiplicity 
of  waveforms.  The  approach  here  is  to  characterize  these  processes  with  a 
mathematical  model  (such  as  a  multichannel  AR  description)  and  to  implement  a 
likelihood  ratio  whose  magnitude  is  sensitive  to  the  difference  between  the  model 
parameters  under  each  hypothesis.  For  passive  detection  applications,  we  view 
the  processes  as  arising  from  internal  physical  mechanisms  which  give  rise  to  the 
emission  of  radiation  which  may,  in  general,  contain  partial  correlation  when 
observed  over  specific  bands.  In  biomedical  applications,  we  might  consider  the 
processing  of  EEG  waves  where  we  may  seek  to  detect  a  weak  brain  potential 
among  other  strong  brain  signals.  For  the  purpose  of  validating  the  theoretical 
results  developed  here,  at  a  later  time,  consideration  will  be  given  to  the 
detection  problem  utilizing  the  data  collected  from  three  co-located  radar 
systems  simultaneously  operating  at  three  distinct  frequencies.  The  resulting 
observation  data  will,  however,  be  treated  as  partially  correlated  across  the  three 
channels.  Performance  evaluations  will  be  determined  in  terms  of  receiver 
operating  characteristics.  Both  analytic  as  well  as  Monte  Carlo  approaches  will 
be  utilized  in  the  analysis. 

A  likelihood  ratio  for  this  detection  problem  with  stationary,  Gaussian 
signal  and  interference  processes  has  been  developed  in  section  V.  In  section  VI, 
implementation  architectures  for  the  likelihood  ratio  are  briefly  presented. 
Section  VII  outlines  the  future  investigation.  At  a  later  time,  we  shall  address  the 
use  of  adaptive  methods  to  consider  the  detection  problem  when  the  statistical 
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processes  are  unknown  and  time  varying.  The  adaptive  procedure  is  also 
approached  by  postulating  an  underlying  parametric  model  for  the  observation 
processes.  This  approach  enables  the  use  of  various  adaptive  algorithms  to 
estimate  the  model  coefficients  and  thus  utilize  these  estimates  to  update  the  filter 
coefficients.  The  adaptive  procedures  are  utilized  to  retain  a  robustness  in  non¬ 
stationary  processes. 
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II.  THE  DETECTION  PROBLEM 


In  the  multichannel  simple  binary  detection  problem,  the  discrete  revived 
baseband  waveforms  are 


Hi: 

i(n)  = 

fi  (n)  +  £  (n)  +  w.  (n) 

n  *  1,2,.. 

•»N 

Hq: 

X  (n)  = 

£  (n)  +  w  (n) 

n  =  1,2,., 

,,N 

(2-1) 

where  a  (n)  is  a  zero  mean,  stationary  Jxl  received  observation  vector  consisting 
of  J  channels  and  &  (n),  £  (n)  and  &  (n)  are  zero  mean,  complex  Gaussian  random 
Jxl  vector  processes  describing  the  signals,  non-white  noise  and  white  noise, 
respectively.  We  will  assume  that  the  white  noise  process  is  uncorrelated  with 
s(n)  and  £  (n)  and  is  furthermore  uncorrelated  with  itself  in  time,  but  not  across 
channels,  so  that 


E 1*  (n)  s'11  (k)]  = 


[0] 


■  ww 


(0) 


n#c 

n=k 


(2-2) 


where  Rww(0)  is  the  JxJ  correlation  matrix  of  win).  The  vector  processes  $(n) 

and  £(n),  however,  contain  an  arbitrary  correlation  in  time  and  between  channels. 
We  will  consider  the  condition  where  $(n),  £(n)  and  atfn)  are  jointly  wide-sense 
stationary  processes.  The  correlation  matrix  for  the  observation  data  expressed 
in  index  ordered  form  [1]  is 

Rx  a  =  E[&i  ,n*i  ,n  1  (2'3) 

where 

=  (lT  (1)  JcT  (2)...£t  (N)]  (2-4a) 

iT(k)  =  (x,(k)  x2  (k)...x,(k)].  (2-4b) 


Under  the  condition  of  stationarity,  Ri  is  a  Hermitian,  positive  semi-definite 
matrix.  Furthermore,  this  matrix  can  be  written  in  block  form  as 
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Under  the  condition  of  stationarity,  R*  is  a  Hermitian,  positive  semi-definite 
matrix.  Furthermore,  this  matrix  can  be  written  in  block  form  as 


R*x(0)  RXX(-1)...RXX(.N+1) 

Rxx(0)  RH^l)...RH»(N.l) 

II 

CQ  8 
06 

Rxx(l)  Rxx(0)...Rxx(-N+2) 

• 

s 

RHxx(-1)  RHxx(0)  -RHxx(N-2) 

• 

3X(N-1)  Rxx(N-2)-  -RicX(0)_ 

Rxx(-N-1)  R“,.(-N+2)...R"x(02_ 

(2-5) 


where 


R*x  0)  =E  (x  (k)  (k-|)]  k  =  1,2,. ...N 

I  =  0,  ±  (N-l)  (2-6) 

and  the  last  expression  in  eq  (2-5)  results  because  Rxx(0  =  RHxx  (-•)•  It  is  noted, 
however,  that  each  block  matrix  of  RM  is  not  Hermitian;  i.e.,  RXx  (I)  *  RHxx  (0 
for  MO.  We  also  note  that  R^  is  block  Toeplitz.  The  superscript  B  denotes  that 
Rji  is  written  in  block  form  where  each  block  as  defined  in  eq  (2-6)  is  a  JxJ 
correlation  matrix  over  the  J  channels. 
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III.  AUTOREGRESSIVE  PROCESS  MODELS 
A.  DEFINITION  OF  THE  AR  PROCESS 

In  this  analysis,  the  multichannel  observation  processes  obtained  under 
hypotheses  Hi  with  i  =  0,  1  are  assumed  to  be  generated  by  multichannel 
autoregressive  processes.  For  the  single  channel  case,  analyses  have  been 
conducted  [26,27]  which  indicate  the  appropriateness  of  such  models  for  radar 
applications.  For  the  multichannel  case,  similar  investigations  remain  as  a 
potential  area  for  future  research. 

The  multichannel  Jxl  vector  process  (nlHj)  with  i  *  0, 1  is  expressed  as 
^  H 

x(nlHj)  =  -  £  A  MiOdHiteOi-k)  +  itfnlHj)  i  =  0,1  (3-1) 

k=l 


where  AMj(klHi)  is  the  kth  JxJ  matrix  coefficient  for  an  AR  process  of  model 

order  Mi.  We  note  that  it  is  expressed  in  terms  of  the  Hermitian  operation  for 
notational  convenience,  but  is  not  treated  here  as  a  Hermitian  matrix.  The  vector 
U.(n)  is  a  Jxl  white  noise  driving  vector  which,  in  general,  has  an  arbitrary 
correlation  across  the  J  channels  so  that 


E  (u  (n)  u VO  = 


r  [0]  1*0 

*-Ruu(0)  1  =  0. 


(3-2) 


Ruu  (0)  is  a  JxJ  covariance  matrix  of  the  vector  process  u(n)  and  may  have  off- 
diagonal  components.  Since  u  (n)  is  uncorrelated  in  time,  but  retains  an  arbitrary 
correlation  across  channels,  then  with  wide-sense  joint  stationarity  of  the  channel 
processes  assumed,  we  can  consider 

U  (n)  =  Cy  (n)  (3-3) 


10 


where  the  J  x  J  matrix  C  is  a  constant  matrix.  This  matrix  gives  rise  to  the 
channel  correlation  on  u  (n).  The  vector  y.  (n)  is  a  Gaussian  white  noise  vector 
unconelated  in  time  and  across  channels  such  that 


E  (s  (n)  * H(n- 


[0] 

Dv 


1*0 


1=0. 


(3-4) 


The  elements  of  the  diagonal  matrix  Dv  are  the  variance  terms  associated  with  the 
white  noise  driving  term  on  each  channel.  And  so,  from  eq  (3-3)  we  can  obtain 
the  zero-lag  correlation  matrix  (assuming  wide-sense  stationarity) 


Ruu(0)  =  E[u(n)uH(n)]  (3-5a) 

=  E[Cy(n)yH(n)CH]  (3-5b) 

=  CDyCH.  (3-5c) 

We  could  assume  unit  variance  on  all  elements  of  Dv  without  loss  of  generality  so 
that  Dv  =  I  and  eq(3-5c)  implies  the  Cholesky  decomposition.  The  significance 
of  this  discussion  is  that  the  correlation  matrix  RUu(0)  is  a  constant  matrix 
associated  with  the  white  noise  driving  term  u  (n).  The  correlation  between  the 
channel  elements  of  u  (n)  gives  rise  to  the  off-diagonal  terms  in  Ruu(0).  It  will  be 
shown  that  this  correlation  causes  the  error  output  vector  £(n)  resulting  from  a 
MMSE  estimation  process  to  retain  some  residual  correlation  across  the  J 
channels.  Since  RUu  (0)  expressed  in  eq  (3-5)  is  Hermitiant,  positive  semi- 
definite,  we  can  perform  an  LDLH  decomposition*  such  that 

Ruu(0) =  LuDuLijII  (3-6) 


where  Lu  is  unit  diagonal  lower  triangular.  Solving  for  Du,  we  obtain 
Du  *  Lu'lRuu  (0)  (Lu"1^ 


(3-7a) 


=  E[V1u(n)uH(n)(Lu*1)H]  (3-7b) 


t  It  is  noted  that  in  general  the  correlation  matrix  Rm,  (0  is  not  Hennitian  for  M0. 
*The  motivation  for  using  the  LDLH  decomposition  is  noted  at  the  end  of  section  IV. 


11 


where 


*  E[£(n)*H  (n)] 

£  (n)  =  W1  u  (n) 


(3-7c) 


(3-8) 


so  that  i  (n)  is  a  J  x  1  vector  containing  uncorrelated  elements.  It  represents  an 
underlying  process  of  the  multichannel  AR  process  which  can  be  viewed  as  a 
"spatially-causal"  white  noise  driving  term.  Since  Lu_1  is  also  lower  triangular 
unit  diagonal,  it  is  invertible  so  that  from  eq  (3-8) 


il  (n)  =  Lu  £  (n). 


(3-9) 


Eq  (3-9)  indicates  that  u  (n) ,  originally  defined  in  eq  (3-3) ,  could  identically  be 
generated  by  the  £  (n)  process  through  the  transformation  matrix  Lu;  i.e.,eq  (3-1) 
can  be  written  in  the  equivalent  form 
3  H 

i(n|Hj)  =  -  £  A  Mj(k|Hj)£(n-k)  +  L^Hj^nlHj)  i  =  0,1  (3-10) 

k=l 

where  Lu  (Hi)  denotes  the  specific  matrix  Lu  under  hypothesis  Hi.  In  section  IV  a 
two  stage  multichannel  prediction  error  filter  is  considered  which  uses  estimates 

of  the  AMj(k|Hi)  coefficients  to  obtain  an  approximation  of  u  (n)  in  the  first  stage 

and  an  estimate  of  Lu~l  to  obtain  an  approximation  of  the  temporally  and  spatially 
uncorrelated  process  £  (n|Hi). 

B.  THE  YULE- WALKER  EQUATION 

The  relationship  between  the  matrix  coefficients  Ajvf(k),  the  covarince 
matrix  of  the  forward  AR  driving  noise  vector  and  the  known  correlation 
matrix  R^  noted  in  eq  (2-3)  can  be  expressed  [1]  as 

-  {(XflM  [0)  -[0)}  (3-11) 


where 
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Am  =  [I  Am(1)  Am(2)  ...  A^M)]. 


(3-12) 


The  matrix  [  Rxx]  is  the  reversed  order  correlation  matrix  of  [R**];  i.e.,  the 
correlation  matrix  obtained  with  the  time  order  of  the  vector  i1>N  from  eq  (2-4) 
reversed.  The  corresponding  equations  for  the  stationary,  backward  AR  process 
is  expressed  as 


Bm^xx^  *  {[0].««[0] 


where 


(3-13) 


(3-14) 


and  [ZblM  is  the  covariance  matrix  of  the  backward  AR  driving  noise  vector. 
Eqs(3-ll)  and  (3-13)  are  the  augmented  forms  of  the  multichannel  Yule-Walker 
equations  and  are  presented  in  more  detail  in  [38]. 
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IV 


DATA  TRANSFORMATION  VIA  LINEAR  PRB1 


The  output  of  a  multichannel  linear  prediction  error  filter  of  order  P  is 
expressed  as 

£(n)  *  i(n)  -  *(n|n-l)  <4-la) 

P  H 

»l(n)+  IAp(k)i(n-k)  (4-lb) 

k-1  -:;- 


»  X  Ap(k)  i(n-k)  (4-lc) 

k=0 

JJ 

where  Ap(k)  k=l,2,...,P  are  the  matrix  coefficients  of  the  linear  predictor, 

Ap(0)  =  I  the  JxJ  identity  matrix,  the  subscript  P  distinguishes  the  matrix 

coefficients  as  belonging  to  a  filter  of  order  P,  and  H  denotes  the  Hermitian 
operation  (i.e.,  the  complex  conjugate  transpose  operation). 

In  Appendix  A,  it  is  shown  that  under  the  condition  that  the  matrix 
coefficients  in  eq  (4-1)  satisfy  the  multichannel  normal  equations, 

E[£(n)£H(n.k)]»[0]  k>0  (4-2) 

and  the  output  vector  process  £(n)  is  a  MMSE  process.  Eq  (4-2)  is  the 
orthogonality  principal  and  indicates  that  the  sequential  outputs  of  the  MMSE 
filter  are  orthogonal  in  time.  The  multichannel  normal  equations  which  are  to  be 
satisfied  to  maintain  this  condition  are  expressed  as 


where 


and 


Ap  (fin)  -  {[£f)p  [0]  ...[0]} 

(4-3a) 

Ap  *  [I  Ap(l)  Ap(2)  Ap(P)] 

(4-3b) 

[ftjsl  -  E[I1N  i1N]  -  EIsn,,  iN.l) 

(4-3c) 

[Ef]P  -  E[E(n)£H(n)]  -  R^fO). 

(4-3d) 
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iftxx]  is  the  reversed  order  known  correlation  matrix  of  eq  (2-3);  i.e.,  the 

correlation  matrix  expressed  with  the  time  order  of  the  vector  in  eq  (2-4) 
reversed.  We  also  note  that  eqs  (4-3a)  and  (3-11)  are  identical  in  form.  This 
equality  implies  that  the  MMSE  estimate  of  the  AR  observation  process  is 
obtained  when  the  prediction  error  filter  coefficients  are  identically  equal  to  the 
AR  process  coefficients.  Figure  4-1  shows  the  synthesis  and  analysis  procedure. 


SYNTHESIS  ANALYSIS 


Figure  4-1 
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With  i(n)  and  u(n)  initially  considered  as  non-random*  processes,  further 
insight  can  be  obtained  by  considering  the  z-transform  of 


M  H 

i(n)  *  -  £  AM(k)x(n*k)  4-  u(n) 
k=l 

so  that 

M  H  c 

X(Z)  =  - 1  AmOO  X(Z)Z*k  +  H(Z). 
k=l 

Bringing  the  summation  to  the  LHS,  we  can  write 


M  TI  | 

X  AM(k)Z'k 

k=0 


X(Z)=U(Z) 


-N44) 

(4-5) 


(4-6) 


where  Am(0)  =  I.  We  now  define  a  filter  representation  for  the  model  process  as 

HM(Z)  = 

so  that 


M  h  u 

X  AmGOZ* 

k=0 


(4-7) 


X(z)  =  Hm(z)  IK. z). 


(4-8) 


Eq  (4-8)  indicates  that  X  (z)  is  the  output  of  the  filter  Hm  (z)  with  input  U  (z). 
Similarly,  the  z-transform  of  eq  (4- lb)  where  &  (n)  and  g  (n)  are  considered  non- 
random  is  expressed  as 


E(Z)  = 


P  H 

X  A?(k)Z-k 
k=0 


X(z) 


(4-9a) 


t  We  initially  consider  non-random  processes  since  the  z-transform  of  a  random  process  is  not 
defined. 
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(4-9b) 

(4-10) 

(4-11) 


When  P£M,  consider  the  case 


aH  Aw(k)  k<M 
AP  = 

l[0]  k>M. 


We  then  have 

HF  (Z)  -  Hm(Z). 


(4-12) 

(4-13) 


At  this  point,  we  can  now  consider  the  input  and  outputs  of  these  filters  to  be 
random.  Using  eqs  (4-8)  and  (4-13)  in  eq  (4-9b) 


E  (z)  =  Hm(z)  Hm  (z)  U(2)  =  u  (z) 


(4-14) 


In  the  time  domain,  eq  (4-14)  is  equivalent  to 

£(n)  =  u(n).  (4-15) 

And  so,  under  the  condition  that  the  prediction  error  coefficients  are 
identical  to  the  coefficients  of  the  AR  model  process  and  under  the  assumption 
that  the  AR  process  is  the  exact  model  of  the  observed  process,  the  prediction 
error  filter  output  g  (n)  is  a  white  noise  vector  equivalent  to  the  AR  model  white 
noise  driving  vector.  However,  it  must  be  emphasized  that  the  use  of  an  AR 
process  with  a  white  noise  driving  function  is  usually  an  approximate 
representation;  i.e.^t  is  not  used  to  describe  the  underlying  physical  mechanisms 
which  give  rise  to  the  random  processes.  Rather,  it  is  a  representation  which  has 
a  system  transfer  function  given  by  eq  (4-7).  We  must  therefore  make  a 
distinction  between  the  model  of  the  processes  (synthesis)  and  the  estimation 
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process  (analysis)[l).  In  general,  the  output  g  (n)  of  the  linear  predictor  is  not  a 
white  noise  vector  output  due  to  the  approximate  representation  of  physical 
processes  by  an  AR  model.  It  is  also  due  to  the  fact  that  we  often  do  not  have  'a 
priori'  knowledge  regarding  the  values  of  the  coefficients  of  this  approximate 
model.  As  a  result,  we  must  estimate  these  coefficients  from  the  observation  data 
as  we  obtain  it.  With  a  limited  amount  of  data,  the  filter  coefficients  are  only 
estimates  of  the  AR  process  coefficients. 

For  stationary  processes,  these  coefficients  could  be  determined  through 
estimates  of  the  multichannel  correlation  matrix  lag  values  and  the  Levinson- 
Wiggins-Robinson  algorithm  [3].  Although,  other  methods  proposed  by  Strand- 
Nuttall  and  Morf- Vieira  have  been  developed  [24,25,30,31]  with  improved 
performance  with  limited  data.  For  non-stationary  processes,  adaptive  schemes 
must  be  considered.  We  will  address  this  topic  in  a  subsequent  report. 

At  this  point,  we  note  that  eq  (4-15)  resulted  from  the  analysis 
procedurefvia  a  linear  prediction  error  filter  with  coefficients  given  by  eq(4-12)] 
of  the  process  synthesized  by  eq  (4-4).  If  u  (n)  is  assumed  to  be  uncorrelated 
across  channels,  the  resulting  £  (n)  is  also  uncorrelated  in  time  and  space  (i.e., 
channels).  In  general,  as  noted  in  the  previous  section,  u.  (n)  may  possess 
arbitrary  correlation  between  the  J  channel  elements.  Therefore,  the  vector  £  (n) 
will  retain  a  residual  correlation  over  the  channels  due  to  the  spatial  correlation 
of  u  (n). 

Since  the  matrix  Ree(0)  =  E  [£  (n)  £H  (n)]  is  Hermitiant ,  and  positive  semi- 
definite,  we  can  perform  an  LDLH  decomposition*  such  that 

ReeW^LyDyLyH  (4-16) 

Solving  for  Dy 

Dy  =  Lyl  Re^OXLy^H  (4-17a) 

=Ly'E  [£(n)  £»(„))  (Ly  1  )H  (4-17b) 

t  It  is  noted  that  in  general  the  correlation  matrix  R^e  (I)  is  not  Hermitian  for  MO. 

*  Other  decompositions  could  have  been  used  such  as  Cholesky  or  unitary(51],  however,  the  motivation  for  the 
LDLh  decomposition  is  noted  at  the  end  of  this  section. 
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=E  [Ly'l£(n)  £H(n)(Ly-l)H]  (4- 17c) 

=H  (n)  ^H(n)]  (4-1 7d) 

where 

X  (n)  =  Ly*1  £  (n)  (4-18) 

so  that  the  vector  x  (n)  contains  uncorrelated  elements.  Also, 

E  [  £  (n)  £»  (n )]  =  [0]  n*n\  (4-19) 


Then,  using  eq  (4-18)  to  solve  for  £(n)  and  substituting  this  result  in  (4-19),  we 
obtain 

E  [Ly  x  (n)  (n')  L^H]  =  [0]  n^  n  (4-20) 

so  that 

Ly E  [y (n)  &  (n))  LyH  »  [0]  n *  n.  (4-21) 

Finally, 

E[^(n)^H(n)]*[0]  n*n\  (4-22) 


Eq  (4-22)  implies  that  x  (n)  retains  its  temporal  decorrelation  while  eq(4-17d) 
denotes  its  spatial  decorrelation.  When 

H  fAM(k) 

AP*  m  -  (4-23) 


l[0]  k>M 

the  output  of  the  first  filter  stage  converges  toward  u(n),  so  that  eq  (4-18) 
becomes  (noting  that  Rec  (0)  »  RUu  (0)  and  the  uniqueness  of  the  LDLH 

decomposition) 


X  (n)  =  Lg-1  £  (n)  *  (n) 

*  Lu‘^Lu  £(n) 
*£(n) 


(4-24a) 

(4-24b) 

(4-24c) 


where  eq  (3-9)  was  used  to  obtain  eq  (4-24b).  Thus,  as  the  filter  coefficients 


converge  to  the  coefficients  of  the  AR  process,  x  (n)  approximates  the  spatially 
and  temporally  whitened  process  z(n).  In  addition,  y(n)  has  been  obtained  through 
a  causal  and  causally  invertible  transformation  of  the  original  observation  process 
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&  (n);  ie.  the  input  vector  *  (n)  could  be  recovered  through  an  inverse  filter 
operation  on  y  (n).  This  "information  preserving"  feature  of  ^  (n)  justifies  its  use 
in  a  likelihood  ratio  test. 

A  procedure  that  could  be  used  to  obtain  the  estimate  Lu*1  consists  of 
estimating  the  correlation  lag  values  of  Ree(O)  using  time  samples  of  £  (n).  An 
LU  decomposition  of  this  matrix  would  provide  Ly.  The  inverse  matrix  Ly~l 

would  then  be  the  required  estimate  t-u'1-  An  alternate  approach  based  on  a 
Gram-Schmidt  procedure  will  be  investigated  using  correlation  lag  estimates  of 
Ree(O).  The  motivation  for  considering  the  LDLH  decomposition  is  based  on  the 

anticipation  of  utilizing  a  single  stage  recursive  procedure  to  obtain  the  filter 
coefficients  required  to  estimate  ^n)  {see  section  VI  and  [36]}. 
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V.  MULTICHANNEL  LIKELIHOOD  RATIO 
A.  DERIVATION 

In  this  section,  we  develop  a  multichannel  likelihood  ratio  for  the  detection 
problem  of  eq.(2-l);ie. 

Hj:  *(n)  « s(n)  +  £(n)  +  ^(n)  n  *  1,2, ...,N  (5-la) 

Hq:  i(n)  *  £(n)  +  ^(n)  n  =  1,2,...,N  (5-lb) 

where  each  of  the  complex  vectors  are  Gaussian,  J  x  1  vectors  and  &(n)  is  a 
baseband  observation  vector.  The  methodology  derived  here  stems  from  the 
considerations  presented  in  [51]  for  real  processes.  Under  hypothesis  Hj,  the 

multivariate  joint  Gaussian  density  can  be  written  as  the  product  of  conditional 
densities  so  that 

N 

Px<Sl,nlHi)  -  p[iO)IH;)]  II  p[iO>lii,n-l.Hil  i=0>>  (5‘2) 

n=2 

where 

J,n>[iT(l)AT(2)...ST(n)J  (5-3) 

*T(k)  =  [Xl(k)  x2(k)...xj(k)]  (5-4) 

and 

ll,l=i(l)  (5~5) 

and  all  the  conditional  densities  are  Gaussian.  The  mean  of  the  multivariate 
conditional  density  p[i(n)|ii>n_i,Hj]  is  &(n|n-l,Hi);  ie.  the  linear  MMSE  predictor 
of  i(n)  using  past  data  n.j  and  assuming  Hj  is  true.  The  JxJ  covariance  matrix 

of  this  density  function  is  the  conditional  covariance  matrix  Kx(n|n-l,Hj)  such 
that 

Kx(n|n-l,Hj)  =  E{(s(n)  -  i(n|n-l,Hj)][s(n)  -  i(n|n-l,Hj)]H}  (5-6a) 

-  EtonlHjXfi^nlHi)]  n  -  1,2,...,N 

i  -  0,1  (5-6b) 
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where 

££nlHj)  =  2L(n)  -  £(nln-l  tH[)  i  =  0,1  (5-7) 

is  the  zero-mean  MMSE  vector.  Assuming  wide-sense  joint  stationarity  on  the 
bandpass  processes  ,  the  conditional  density  functions  can  be  expressed  in  terms 

of  the  quadratic  form*  such  that 
p(&(n)l&]>n.j  Jiil  = 

-  — j— ~  1 - exp{[i(n)  -  i(nln-l,Hi)][Kx(nln-l,Hj)]'1[i(n)  -  i(nln-lJfij)]H} 

(n)JIKx(nln-l,Hj)l 


Using  eq(5-7)  in  (5-8),  we  obtain 


i  =  0,1.  (5-8) 


p[*(n)lxj  j.Hj]  = 

=  —i - - - exp  { -  £H(nlH:)[Kx(nln- 1  Mi )]'  ^(nlH;) }  i  =  0, 1  (5-9) 

(7i)JIKx(nln- 1  ,H  j)l  1  X  1  1 

where 

pUCUlHj]  =  plfiCUlHj]  = 

-WJ|KxOIO.Hi)ieXp(-eH(1'Hi)1K-(llHi)1'l6(1'Hi)l  i‘<U  (5'10) 

and 

A(ll0,Hj)  =  0  i  =  0,1.  (5-11) 

We  can  now  express  the  likelihood  ratio  for  the  multivariate  joint  Gaussian 
density  functions  as 


*  The  condition  of  wide-sense  stationarity  provides  specific  relationships  between  the  auto-  and  cross-correlation 
functions  of  the  in-phase  and  quadrature  components.  These  relationships  enable  the  multivariate  Gaussian  density 
function  to  be  expressed  directly  in  terms  of  the  complex  random  vectors  as  given  by  eq(5-8). 
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(5-1 2a) 


Px&l,NIHl) 

HiiHo'Px(Ai,n'H0) 

N 

(p[x(l)IHi]  n  p[x(n)lx.i,n-l»Hl]) 

= - ^ -  (5 -12b) 

(p[x(1)IHq]  n  pfi(n)lii  n-l»^Ol) 

n=2 


where  the  last  equality  results  from  eq(5-2).  Substituting  eqs(5-9)  and  (5-10)  into 
(5- 12b)  and  taking  the  natural  logarithm,  we  have 


InA 


H,  ,H 


l’n0 


n  IKx(nln- 1  ,H0)I  exp  [  -£H(nlH x  )[Kx(nln- 1  ,H  j )]” 1 1£  { nIH  { ) ) 
n=l 

_N  H  -1 
FI  IKx(nln- 1  ,H  j )l  exp [  -£H(nlH0)[Kx(nln- 1  ,H0)]  £(nlH0) } 

_n=l 


N 

=  I 

n=l 


IKx(nln-l,H0)l  , 

lnIKx(nin-l,Hi)l +  £H("'Ho)[Kx(nln- 1  ,H0)]'  *£  ( nlH0) 


(5-13a) 


-  £H(nlH !  )[Kx(nln- 1  ,H  j )]'  *£(nlH  j )]  _ 


(5-1 3b) 

Eq  (5- 13b)  can  be  simpified  further  by  a  diagonalization  of  the  conditional 
covariance  matrices.  Since  these  matrices  are  Hermitian,  we  can  perform  their 
LDLh  decomposition  such  that 

Kx(nln-t,Hj)  =  Ly(nlHj)  D^nlHj)  L^(nlHj)  i  =  0,1  (5-14) 


where  the  matrix  D-y(nlHj)  i=0,l  is  a  diagonal  matrix  with  elements  dX(nlHj).  In 

the  discussion  below,  we  note  that  Kx(nln-l,Hj)  is  a  deterministic  quantity[59]. 

We  therefore  consider  this  matrix  as  well  as  those  on  the  RHS  of  eq(5-14)  as  non- 
random.  Solving  for  D^(nlHj),  we  have  from  eqs(5-14)  and  (5-6b) 
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D^nlHj)  =  (n|Hi)Kx(n|n- 1  ,Hj)[L^  (nlHj)]" 1  i  =  0,1 

«  L^OilHj)  E[£(n|Hi)(£H(n|Hi)]  [L^(n|Hi)]H 

=  eJ 


Ly\n\Hi)  £(n|Hi)£H(n|Hi)  [L;1(n|Hi)]H 


=  EWnlHi^nlHi)] 


i  =  0,1 


(5- 15a) 

(5- 15b) 

(5- 15c) 
(5-15d) 


where 


-1 


(5-17) 


(5-1 8a) 


^<n|Hj)  =  Uy  (nlH^nlHi)  i  =  0,1.  (5-16) 

Eq(5-15d)  implies  that  2(n|Hj)  contains  uncorTelated  elements  across  the  channels. 

From  the  orthogonality  condition,  we  have 

E[£(n|Hj)£H<n'|Hi)  =  [0]  n  in'  i  =  0,1. 

Solving  for£(n|Hj)  in  eq(5-16),  eq(5-17)  becomes 

E{L7(n|Hi)  jdilHj^di’lHi)  [L^(n'|Hj)]H}  =  [0] 

n  ±  n  i  =  0,1 

so  that 

LydilHj)  EttfnlHjJjVlHj))  [L^(n'|Hj)]H  -  [0] 

n  ^  n’  i  =  0,1. 

Since  the  matrices  Ly(n|Hj)  are,  in  general  non-zero,  it  follows  that 
E(y(n|Hj)yH(n'|Hi)]  =  [0]  n  ±  n'  i  =  0,l. 

In  summary, 

E^<"IHi)lH(”'IHi)]  =  {(D0f|Hi)nn:", 

From  eq(5-14),  we  now  have 

K^nln-l.Hj)  =  [L^nlH^l^D'VnlH^L’^nlHi)  i=0,l 


i  =  0,1. 


(5-1 8b) 
(5-19) 

(5-20) 

(5-21) 
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and 


|Kx(n|n-l,Hi)|  =  |Ly(n|Hi)||D7(n|Hi)||L“(n|Hi)|  i=0,l 

=  P^nlHi)! 

J  J  y 

=  II  d^nlHj)  =  riS(nlHi)  i=0,l 

H  j=i 


where 

and 


|L-y(nJHj)|  = 


=  1 


i  =  0,1 


dT  (nlHj)  =  E[|Yj(n!Hj)r]  =  o^j(n|Hi)  i  =  0, 1 . 


(5-22a) 

(5-22b) 

(5-22c) 


(5-23) 

(5-24) 


-1 


The  quantity  a^,j(n|Hj)  denotes  the  variance  of  Yj(n|Hj).  Using  eq(5-21),  the 
quadratic  terms  in  eq.(5-13b)  become 

£H(n|Hj)[Kx(n|n-l,Hi)]'*£(n|Hj)  = 

=  £H(n|Hi)[I^(n|Hi)]-lD‘1(ii 

=  [L^'(n|Hj)£(n|Hj)]H  Dy‘(n| 

=  ^(nlHjJD^  (n|Hj)j(n|Hj) 

f 

J  lYj(n|Hj)|2 

=  2 
j=l  o^(n|Hj) 

where  the  last  two  equations  result  from  eqs(5-16)  and  (5-24),  respectively. 
Using  eqs(5-22c)  and  (5-25d)  in  eq(5-13b),  we  obtain 


i  =  0,1 


i=0,l 

(5-25a) 

i=0,l 

(5-25b) 

i  =  0,1 

(5-25c) 

(5-25d) 
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J  N 

lnAH„Ho  -  1  X 
j=l  n=l 

Eq(5-21)  is  the  multichannel  likelihood  ratio  for  wide-sense,  joint 
stationary,  Gaussian  processes.  It  represents  a  generalization  of  the  single  channel 
likelihood  ratio  reported  in  [6,7,23,39]  for  Gaussian  processes.  In  section  IV,  we 

discussed  a  two  stage  filtering  method  using  multichannel  linear  prediction 
filtering  to  obtain  y(n|Hj).  For  a  known  correlation  function  under  each  of  the 

two  hypotheses,  two  sets  of  filter  coefficients  and  error  variances  can  be  obtained 
exactly,  through  the  multichannel  Yule-Walker  equations.  In  this  case,  y(n|Hj) 

will  be  a  MMSE  output  of  the  filter  for  which  the  hypothesis  Hj  is  true.  The  other 

filter  output,  however,  will  increase  in  terms  of  its  average  output  magnitude. 
Eq(5-21)  indicates  that  when  Hq  is  true,  the  last  term  will  contribute  a  larger 

value  in  a  negative  sense,  causing  the  likelihood  ratio  to  decrease.  Alternatively, 
when  Hi  is  true,  the  second  term  increases  in  a  positive  sense  so  that  the 

likelihood  ratio  increases.  For  the  unknown  correlation  case,  the  filter 
coefficients  and  error  variances  must  be  estimated.  In  this  case,  eq(5-21)  becomes 
a  generalized  likelihood  ratio  and  is  therefore  suboptimal.  From  one  set  of 
observation  data,  we  must  estimate  the  parameters  (ie.  the  filter  coefficients  and 
error  variances)  for  each  filter  assuming  the  appropriate  hypothesis  is  true.  In  the 
practical  implementation  for  this  case,  we  must  assume  that  the  MMSE  filters 
under  each  hypothesis  have  distinct  orders  so  that  the  likelihood  ratio  will,  in 
general,  have  a  value  other  than  zero.  This  assumption  is  justified,  for  example, 

when  characterizing  the  observation  data  as  autoregressive  processes  where  the 
order  of  the  process  under  hypothesis  Hj  (signal  present)  is  larger  than  that  under 

Hq  (signal  absent).  For  single  channel  processes,  these  considerations  have  been 
treated  in  [6,8]. 


In 


j  (n|H0)  lYj(n|H0)|2  l7j(  n|  H ; 


J 

Oj”(n|H|) 


of(n|H0) 


(5-21) 
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B.  THE  SPECIAL  CASE  OF  INCOHERENT  INTEGRATION 

In  [59],  several  limiting  forms  of  the  single  channel  likelihood  ratio  are 
discussed.  They  include  (1)  the  detection  of  a  deterministic  signal  in  additive 
white  Gaussian  noise,  (2)  the  detection  of  a  non-white  Gaussian  signal  process  in 
additive  white  Gaussian  noise,  (3)  the  moving  target  detection  (MTD)  algorithm 
for  the  detection  of  a  deterministic  signal  in  non-white  Gaussian  noise  [60,  61] 
and  (4)  the  detection  of  a  non-white  Gaussian  signal  process  in  additive  non-white 
Gaussian  noise  [62].  It  is  also  shown  [59]  that  this  likelihood  ratio  contains  the 
algorithm  developed  [63]  for  the  detection  of  a  deterministic  signal  in  additive 
non-white  Gaussian  noise  of  unknown  correlation  statistics. 

One  limiting  form  of  this  likelihood  ratio  which  has  not  been  noted  is  the 
case  of  uncorrelated  signal  processes  in  additive  white  noise.  In  this  case,  the 
signal  is  a  white  noise  process.  Since  past  values  of  such  a  process  are 
uncorrelated  with  present  and  future  values,  the  coefficients  for  the  MMSE 
estimation  are  zeroes  in  the  single  channel  case.  Thus,  the  error  signals  in  this 
case  are  identical  to  the  observation  processes.  For  J=1  and  known  variance 
terms,  we  have 


O2(H0)  =  <4 


(5-22a) 


and 


a2(H !)  =  a^  +  a2 


w 


(5-22b) 


2  2 

where  o  and  o  are  the  known  variances  associated  with  the  signal  and  white 


noise,  respectively.  Absorbing  the  constant  term  into  the  threshold,  eq(5-21) 
becomes 

N 

I  I  TL  I  I  I  II  I  &  I  II  II  I 

(5-23a) 


"■Ah, .h„-  I J 

n=l 


w 


CTs  +  a  w 
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(5-23b) 


2 

°s  N 

=  ~2 — 2 - 1~  Zlx(n)'  • 

ow(as  +  ow)  n=l 

Eq(5-23b)  is  the  likelihood  ratio  for  complex,  uncorrelated  signal 
processes  in  additive  white  noise.  Its  form  is  noted  in  [12, Part  I, Chap.  2]  for  the 
case  of  real  processes.  Eq(5-23b)  is  also  informative  from  another  point  of  view. 
We  recognize  that  it  has  the  form  of  an  incoherent  integrator.  In  section  VI,  we 
discuss  the  implementation  of  eq(5-21)  with  a  dual  filter  strategy,  where  each 
filter  is  designed  for  the  hypothesis  Hj  i=0.1.  Eq(5-23b)  indicates  that  if  the 

coefficients  of  these  filters  are  set  to  zero,  the  likelihood  implementation  becomes 
a  simple  incoherent  integrator.  More  importantly,  this  consideration  reveals  the 
useful  role  of  the  filtering  process;  ie.,  when  the  filters  are  used  as  prediction 
error  filters,  they  provide  the  means  to  obtain  coherent  integration.  This 
integration  gain  is  achieved  when  the  signal  process  contains  pulse-to-pulse 
amplitude  and  phase  correlation.  The  detection  performance  for  various  signal 
correlation  levels  will  be  described  in  forthcoming  reports. 
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VI.  LIKELIHOOD  RATIO  IMPLEMENTATION  SCHEME 

A.  SYSTEM  ARCHITECTURE 

A  block  diagram  of  a  system  that  approximates  the  likelihood  ratio  of  eq 
(5-22)  is  shown  in  Fig  6-1.  It  is  recognized  as  the  multichannel  extension  of  the 
implementations  reported  in  [8,23].  The  specific  choice  of  the  prediction  error 
filter  structure  will  depend  on  the  assumed  underlying  model  of  the  observation 
process  [14].  A  forward  prediction  error  filter  using  a  tapped  delay  line 
architecture  is  shown  in  Figs.  6-2  and  6-3.  The  lattice  structure  which  utilizes 
both  forward  and  backward  coefficients  is  shown  in  Figs.  6-4  and  6-5. 

We  note  that  the  architecture  in  Fig  6-1  utilizes  two  prediction  error  filters 
(PEF)  implemented  in  parallel  with  each  designed  to  be  an  optimal  estimator 
under  the  given  hypothesis.  As  noted  in  section  V.A,  a  linear  filter  of  a  higher 
order  may  be  used  on  the  filter  selected  under  Hi  as  compared  to  that  for  Ho  if 
the  underlying  process  is  assumed  to  be  an  AR  process  of  a  higher  order  with  the 
signal  present  [8,23].  Under  this  condition,  when  hypothesis  Ho  is  true  (i.e.,  no 
signal  present),  both  filters  provide  a  MMSE  output  when  the  optimum  filter 
coefficients  are  determined.  When  Hj  is  true,  the  filter  designed  for  Hi  adjusts  to 
the  new  process.  However,  the  filter  for  Ho  produces  a  much  larger  error  output 
since  this  lower  order  filter  now  operates  on  a  higher  order  AR  process  and 
therefore  cannot  adapt  to  the  underlying  process  coefficients.  As  a  result,  the 
second  summation  in  eq  (5-22)  increases.  Since  this  term  provides  a  positive 
contribution  to  the  likelihood  ratio,  it  is  the  mechanism  which  raises  the 
likelihood  above  a  predetermined  threshold  under  Hi. 

B.  FILTER  COEFFICIENT  DETERMINATION 

The  determination  of  the  filter  coefficients  can  be  realized  via  several 
approaches.  For  stationary  processes,  the  most  straight  forward  would  be  a 
solution  of  the  Ap(k)  coefficients  via  the  Levinson-Wiggins-Robinson  (LWR) 

method;  howevei,  the  Strand-Nuttall  and  the  Vieira-Morf  methods  may  also  be 
considered  [24,25,29,30,31].  For  the  LWR  method,  estimates  of  the  correlation 
lag  values  as  well  as  the  prediction  error  variances  are  obtained  by  recursive 
update.  Since  this  method  utilizes  correlation  lag  value  estimates,  it  is  anticipated 
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to  yield  less  accurate  estimates  of  the  coefficients  than  the  Strand-Nuttall  or 
Vieira-Morf  approaches  for  the  same  reasons  as  noted  by  Burg  [35]  in  the  single 
channel  case;  ie,  the  unbiased  estimate  of  the  correlation  matrix  may  not  be 
positive  semi-definite  and  therefore  physically  unrealizable.  On  the  other  hand, 
the  biased  estimate  may  yield  inaccurate  estimates,  especially  for  limited  data. 
The  Strand-Nuttali  method  (the  multichannel  generalization  of  the  Burg 
algorithm),  however,  estimates  the  filter  coefficients  directly  from  the  data  thus 
bypassing  the  requirement  to  estimate  the  correlation  matrix.  Actually,  the 
reflection  coefficients  for  the  single  chacael  or  multi-channel  cases  are 
determined  by  this  algorithm  and  applied  to  the  lattice  filter  structures  shown  in 
Figures  (6-4)  and  (6-5),  respectively.  From  these  values,  autoregressive 
coefficient  matrices  AH(k)  [as  well  as  the  backward  coefficients]  can  be 
determined. 

The  additional  matrix  coefficient  Ly1  discussed  in  section  IV  must  also  be 
determined  in  order  to  use  the  simplified  form  of  eq(5-21).  This  coefficient 
completes  the  spatial  whitening  of  the  £  (n)  process  via  eq(4-18).  As  noted  in 
section  IV,  the  elements  of  this  matrix  could  be  obtained  by  performing  estimates 
of  the  correlation  matrix  Ree(O)  =  E  [£  (n)  £H  (n)]  using  the  output  £  (n)  from  the 

first  stage  of  processing.  An  LU  decomposition  of  this  matrix  would  yield  the 
lower  triangular  unit  diagonal  Ly  matrix  which  is  an  estimate  of  L^.  In  [36],  a 

recursive  procedure  is  being  considered  which  computes  the  filter  coefficients  for 
the  structure  shown  in  Figure  6-3.  This  single  stage  filter  is  equivalent  to  that 
shown  in  Figure  6-2;  however,  these  coefficients  are  obtained  via  a  single 
recursive  procedure. 

For  the  LWR  method,  the  unbiased  cross-correlation  function  at  lag  I  between 
channel  i  and  j  at  the  Nth  data  sample  is  calculated  as 

N  i  N 

l^ij  (I)  =  7  £  *i(k)x  j  (k-l)  ij  -  1,2 . J 

k=l+l  1 


I  =  0,1,2 . P  (6-1) 


where  P  is  the  filter  order  and  N  is  the  total  length  of  the  time  series  with  N»P. 
Since  we  are  considering  jointly  stationary  processes  when  using  this  approach, 
we  also  have 
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(6-2) 


N  N  * 

<£<o  -  [^(-i)]  • 

Eqs  (6-1)  and  (6-2)  enable  us  to  fill  the  entire  correlation  matrix  of  eq  (2-3)  with 
updated  estimates.  The  expression  in  eq  (6-1)  can  be  made  computationally 
efficient  by  re-expressing  ii  in  terms  of  the  recursive  form 


/i\  i 

Tij  (')  =  TTT" 


X.(N)  x*  (N-l)  +  -j-  ^  Xj(k)  x*  (k-l) 

J  m  k-ui  ‘ 

x  ,(N)  x^  (N-l)  rN+l  1  i  ^ 

FTi  +  lT5T^j  nTT  ^  i(K,xj  (  ' 

k-l+1 


x  ,(N)  x*  (N-l)  r  N-l-1  -1  AN., 

= — M - +  brrJ  Fij  0) 

Likewise,  the  biased*  estimate  of  rij  (I)  is  expressed  as 


rt  ©  =  TJ-  X  xi(k)x*  (k'l) 

k-  1+1 


(6-3a) 


(6-3b) 


i,j  —  1,2,...,J 
I  =  0,1,...,P.  (6-3c) 


(6-3d) 


x  j(N)  x*  (N-l)  r  N-l  1  AN.,  A  i  i  _  ,  ,  , 

=  - jr* -  +  L — fj —  J  rij  v)  ‘•J  l’z . 3 

1  =  0,1 . P.  (6-3e) 

Eqs  (6-3c)  and  (6-3e)  enable  the  current  estimate  of  rij  (I)  to  be  made  in  terms  of 
past  estimates  rij  (I)  without  a  complete  recalculation  using  all  prior  data. 


*  The  biased  estimate  of  rjj(l)  is  most  often  used  since  it  ensures  positive  semi -definiteness  of  the  correlation 
matm[25,35]. 
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C.  ERROR  VARIANCE  DETERMINATION 


The  unbiased  estimate  of  the  prediction  error  variance  for  the  jth  channel 
under  Hi  for  i  =  0,1,  at  the  Nth  data  sample  could  be  calculated  using  the  sample 
variance;  ie. 

N  2 

(N|Hj)  =—  X  I  Vj  (k|Hj)  -  y.  (N|Hj)|  i  -  0,1  (6-4a) 

where 

—  1  N 

7j  (N|Hj)  =  —  2,  Tj  OlHj)  i  =  0,1.  (6-4b) 

1-1 


A  recursive  form  for  <jj2  (N|Hi)  can  be  expressed  as  (see  Appendix  C) 


Oj  (N-1IH)  +-L.  lY/NIHj 


)l 


-jq-  [vj(N|H)yj*(N-l|H,)+7* (N|H,)  (N-I|H,)] 


+  -7T  iTjCN-HHi)!2. 

N  J  (64c) 

A  very  informative  discussion  regarding  the  effect  of  error  variance 
estimates  on  the  single  channel  likelihood  ratio  is  presented  in  [14]. 
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VII  PROPOSED  INVESTIGATION 


A.  STATIONARY  GAUSSIAN  PROCESS  SYNTHESIS 

In  this  part  of  the  report,  a  future  investigation  is  proposed.  We  will 
assume  the  presence  of  stationary,  Gaussian  processes  for  both  signal  and  noise. 
The  analysis  will  be  further  divided  into  two  areas:  (1)  signal  in  additive  white 
noise  and  (2)  signal  in  additive  non-white  (clutter)  plus  white  noise.  The 
principal  objective  will  be  the  determination  of  receiver  operating  characteristic 
(ROC)  plots  in  terms  of  probability  of  detection  (Pd),  probability  of  false-alarm 
(Pfa),  the  signal-to-noise  (S/N)  and  signal-to-clutter  (S/C)  ratios  associated  with 
multichannel  processes.  The  analysis  will  also  consider  the  number  of  channels 
and  the  number  of  pulses  as  parameters.  Performance  comparisons  will  be  made 
between  the  single  channel  case  and  the  multichannel  case  as  well  as  to  an  analytic 
evaluation  described  later  in  this  section. 

The  synthesis  of  the  random  observation  processes  under  Hj  for  i=0,  1  is 
now  addressed.  These  processes  will  be  utilized  in  the  performance  of  a  Monte- 
Carlo  analyses  for  the  determination  of  the  ROC  plots.  Two  methods  will  be 
considered. 

METHOD,! 

In  the  first  approach  we  can  characterize  the  observation  processes  as 
multichannel  AR  processes  under  each  hypothesis.  We  then  have 
^  H 

A(nJHj)  =  -  X  A  Mj(k|Hj)*(n-k)  +  u(n|Hj)  i  =  0,1  (7-1) 

k=l 

H 

where  Mi ,  AM^klHj)  and  u  (n|Hj)  denote  the  model  order,  the  matrix  coefficients 

and  the  white  noise  driving  term  under  each  hypothesis,  respectively.  In  the 
single  channel  radar  problem,  analyses  have  been  conducted  to  model  radar 
clutter  with  AR  processes  of  relatively  low  order  [20,21].  It  has  also  been  noted 
[8,19,22]  that  the  sum  of  two  AR  processes  yields  an  autoregressive  moving 
average  (ARMA)  process  which  in  turn  can  be  modelled  as  a  higher-order  AR 
process.  Therefore,  assuming  that  the  signal  and  clutter  noise  are  eacl > 
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characterized  by  an  AR  process,  we  would  expect  that  Mi  >  Mo*  For  single 
channel  processes,  values  of  Mo  =  1  or  2  and  Mi  =  4  or  5  have  been  n  ported 
[81  for  radar  applications.  The  extension  of  this  work  to  the  multichannel  case 
remains  an  open  area  of  research.  Eq  (7-1)  could  be  utilized  to  generate  the 
processes  under  each  hypothesis  using  predetermined  values  for  the  coefficients. 
This  approach  would  be  useful  in  the  diagnostics  of  the  filtering  scheme;  i.e.,  one 
could  validate  that  the  filter  coefficients  converge  to  the  known  preassigned 
model  coefficients  as  well  as  assess  the  convergence  rate  and  final  error  variance. 
This  approach,  however,  does  now  allow  control  over  the  variations  of  the 
signal-to-noise  (S/N)  and/or  signal-to-clutter  (S/C)  ratios  for  parametric 
performance  evaluations. 

METHOD  2 


In  the  second  approach,  we  consider  the  multichannel  extension  of  the 
method  suggested  in  [23].  A  complete  description  of  this  approach  is  presented  in 
[74].  This  is  shown  in  Figure  7-1  where  we  generate  the  Jxl  vectors  s(n)  and  £(n) 
as  distinct  multichannel  AR  processes. 


H 


Us(n)_ 


SIGNAL 

Hs(z) 


Ho{ 


11  in) 


w(n) 


CLUTTER 

Hc(z) 


s(n) 


c(n) 


*(nlH])  =  $(n)  +  c(n)  +  yy(n) 
2L(nH^)=£(n)  +  iy(n) 


Figure  7-1 

The  vectors  us(n),  uc(n)  and  w(n)  are  zero-mean,  Gaussian  white  noise 
vectors  uncorrelated  in  time  but  have  an  arbitrary  correlation  across  channels  as 
described  in  sections  II  and  III.  In  this  case, 
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(7-2) 


H[us(n)u^n-I)]=  |s 
E  [uc(n)ii”(n-l)]=  | 
E  [w  (n)  u^(n-l)3=  | 


[0] 

1*0 

Ruu(0) 

1=0 

[0] 

1*0 

cRuu(0) 

1=0 

[0] 

1*0 

Rww(0) 

1=0 

(7-3) 


(7-4) 


where  SRUU(0)  and  cRUu(0)  are  the  J  x  J  correlation  matrices  of  the  white  noise 
driving  vectors  (n)  and  lie  (n),  respectively-  The  JxJ  matrix  Rww(0)  is  the 

corresponding  matrix  for  the  additive  white  noise  vector  $:(n).  In  general,  they 
have  off-diagonal  components.  The  functions  Hs  (z)  and  Hc  (z)  are  the  filter 
representations  for  the  synthesis  of  the  signal  and  clutter  processes,  respectively. 
For  the  generation  of  signal  AR  processes  we  would  use 

Ms 

S(n)  =  -X  Ag(k)  i(n-k)  +u  s(n)  (7-5) 

k=l 


where  M$  is  the  order  of  the  signal  model  and  AsH(k)  is  the  matrix  coefficient  of 
the  process. 

We  must  now  establish  a  procedure  to  determine  the  matrix  coefficients 
H 

As  (k)  which  when  used  in  eq  (7-5)  will  yield  realizable  values  for  s(n).  The 

cross-correlation  functions  for  the  signal  and  non-white  noise  processes  are  to  be 
selected  using  the  functional  forms  [74]: 


fg(Xgij,  I  -  lgij)  1 1=0 


g  =  S,c 


(7-6a) 
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. 


where  g=s,c  refers  to  signal  and  non-white  noise,  respectively;  pe..  is  the  cross 
correlation  coefficient;  Og..  and  <5gjj  are  the  standard  deviations  for  the  channel  i 
and  j  processes;  Xg.j  is  the  temporal  cross-correlation  coefficient  and  is  a  measure 
of  the  correlation  between  pulses  (relative  to  Ig.j)  across  the  channels  i  and  j  for 
i^j  f74];  for  i=j,  Xg..  is  a  temporal  autocorrelation  coefficient  for  the  pulse-to- 
pulse  correlation  on  channel  i  [23];  0g.j(l)  is  the  phase  of  the  cross-correlation 

function  and  allows  for  the  modeling  of  processes  with  uneven  spectral  shape 

including  Doppler  shifts.  The  functions  f(?ie..,l-l0..)  are  selected  to  appropriately 

*ij 

shape  the  autocorrelation  (i=j)  and  cross-correlation  (i  ^  j)  functions.  The 
parameter  Ig.j  is  the  lag  value  for  which  the  function  f(*)  has  a  peak  value  of 

unity  and  accounts  for  the  fact  that  the  cross-correlation  function  does  not  peak  at 
1=0  as  the  autocorrelation  function  does.  Eq(7-6a)  could  be  modified  using  the 
relation 

lpgij'  =  pgij  exP[-jQgij(°)]-  0 -6b) 

For  the  autocorrelation  function  (i=j),  we  note  that  lpg..l=l,  lg..=0  and 

0g..(0)=0  [see  74]  so  that  eq(7-6)  reduces  to  the  expression 

Rfi(l)  =  agii  fg(^gii’|)exPU(0gii(l)}}  g  =  s.c.  (7-7) 


In  reference  [74],  the  selection  of  the  fg(*)  functions  as  well  as  restrictions 

on  the  correlation  function  values  are  discussed  to  ensure  conditions  such  as  the 
positive  semi-definiteness  of  the  correlation  matrix. 

The  following  procedure  is  then  used  to  generate  the  time  sequenced  values 

1231: 


(1) 

(2) 

(3) 


the  desired  shapes  of  the  autocorrelation  and  cross-correlation  values 
are  obtained  using  the  functional  forms  above. 

the  order  of  the  AR  process  (for  synthesis)  is  selected  based  upon 
requirements  to  fit  the  desired  spectrum, 
s  c 

the  values  of  R -(I)  and  R -(I)  are  used  to  form  the  matrix  elements  of 


Rss  and  Rcc- 
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(4)  the  multichannel  Yule-Walker  equations  are  solved  to  determine  the 
matrix  coefficients  Ag  (k)  for  k  =  1,  2,  i.e.. 


Ag  [Rjj]  =  ([If]"  [0]...[0]) 

where 

Ag  =  (I  Ag  (1)  Ag  (2)...Ag(Mg)] 

(If]g  =  E[UgUg] 

iij  =  (1JT<1)  UT(2)...uT(N>] 
uT(m)  =  [u^m)  u2(m)  ...uj(m)] 


g  =  s,c  (7-8) 

g  -  s,c  (7-9) 

g  =  s,c  (7-10) 

g  =  s,c  (7-11) 

m  =  1,2 . N.  (7-12) 


(5)  the  values  of  AsH(k)  are  now  substituted  into  eq  (7-5)  and  used  in  the 
generation  of  s(n). 

The  vector  2c(nlHi)  i=0,  1  is  thus  obtained  as  described  in  Fig  7-1. 


B.  ANALYTIC  PERFORMANCE  ANALYSIS 


With  the  matrices  R^  and  R^  specified  using  elements  determined  from  eq 
(7-6),  it  may  be  possible  to  proceed  with  an  analytic  determination  of  the 
probability  of  detection  (Pd)  and  probability  of  false  alarm  (Pfa).  For  the  single 
channel  case,  analytic  expressions  for  these  quantities  are  discussed  in  [23].  In  our 
notation 


j  (— 4r-Hexp (- 
1  det  (I-2jt  (Rj+ 


-£0-1] 


2£(Rd+  RS)]*R*] 


and 


J°°Re 


r(--^-)[exp(-ffl)-l] 

l  det  [I-2j4  RjR] 


(7-13) 


(7-14) 
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where  T  is  the  threshold  and  the  single  channel  correlation  matrices  are  expressed 
as 


"d  C'  *'w 

(7-15) 

R=  [Rj1-  (Rd+  Rs ) 1  ] 

(7-16) 

Rs=E(siH] 

(7-17) 

Rc=E[££H] 

(7-18) 

&T=[S(1)  S(2)...S(N)] 

(7-19) 

£t=[C(1)  C(2)...C(N)] 

(7-20) 

Equations  (7-13)  and  (7-14)  will  be  considered  to  determine  if  they  can  be 
generalized  to  the  multichannel  case  using  the  matrices  R^,  Rgc  and  Rww. 
Receiver  operating  characteristic  curves  can  then  be  generated  in  terms  of  Pd 
versus  Pfa  for  the  parameters  noted  in  eq(7-6). 
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C.  MONTE  CARLO  ANALYSIS  (STATIONARY  PROCESSES) 

The  simulation  procedure  discussed  in  VII-A  will  be  used  to  synthesize 
multichannel  observation  processes  under  hypotheses  Ho  and  Hi .  A  sequence  of 
N  data  samples  will  be  generated  for  the  parameters  specified  in  eq(7-6).  The 
threshold  level  T  for  a  given  trial  will  be  established  by  generating  x(nlHo)  using 
ns  samples  where  ns  satisfies  the  condition!  17)  ns  >  10/Pfa.  The  simulation  will 

then  be  rerun  to  generate  N  observation  processes  &(nlHi)  as  inputs  to  the 
likelihood  ratio.  For  stationary  processes,  the  filter  weights  will  be  determined 
as  discussed  in  section  VI  B.  The  processes  £(n)  or  y(n)  will  then  be  computed  in 
terms  of  the  likelihood  ratio  [see  eq  (5-21)  as  well  as  Fig.  6-1).  The  number  of 
threshold  crossings  shall  be  used  to  compute  Pd- 

Since  the  number  of  samples  required  to  establish  a  given  Pfa  varies 
inversely  with  Pfa,  low  Pfa  levels  will  require  long  computer  run  times. 
Therefore,  the  Monte  Carlo  approach  may  have  to  be  limited  to  values  such  as  Pfa 
-  103  with  the  resulting  analyses  used  to  confirm  the  analytic  solutions  of  section 
VII  -  B. 

D.  POSSIBLE  EXTENSIONS 

The  generalization  to  non-Gaussian  processes  can  follow  the  approach 
utilized  in  [4,6].  For  the  single  channel  case,  Metford  formed  a  process 
consisting  of  a  variance  normalized  partial  summation  of  innovations  processes. 
In  [6],  he  proved,  (through  a  lengthy  algebraic  manipulation),  that  this  process 
satisfied  a  central  limit  theorem  approaching  a  Gaussian  distribution  with  a  rate 
of  convergence  of  N*1/2  (where  N  is  the  number  of  pulses).  It  was  shown  that 
the  resulting  likelihood  ratio  was  identical  in  form  to  the  single  channel 
innovations  based  detection  algorithm  (IBDA)  established  for  Gaussian  processes, 
but  with  the  additional  requirement  of  processing  over  N  pulses;  i.e.,  the 
likelihood  ratio  converged  to  the  IBDA  with  a  rate  of  convergence  1/VN  for 
large  N.  The  exact  size  of  N  was  not  investigated  specifically  in  [6],  but  in  [4]  it 
is  indicated  that  these  processes  are  expected  to  be  Gaussian  even  for  "a  small 
number  of  samples." 

In  the  previous  sections,  stationary  observation  processes  were  assumed. 
The  condition  of  stationarity  will  be  maintained  for  the  major  portion  of  the 
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performance  evaluations.  However,  the  extension  to  adaptive  methods  to  update 
the  filter  coefficients  can  be  considered.  These  methods  are  important  when 
considering  non-stationary  processes  with  changing  statistics.  Consideration  will 
be  given  to  methods  such  as  the  least-mean-square  (LMS)  algorithm,  recursive- 
least-squares  (RLS)  [70]  and  adaptive  lattice  filter  methods. 

The  potential  for  detection  performance  improvement  when  utilizing  ’a 
priori'  information  to  predetermine  the  filter  coefficients  under  hypothesis  Ho, 
could  be  considered.  In  this  case,  reference  data  collected  on  the  non- white  noise 
processes  could  be  utilized  to  preset  the  filter  weights.  In  the  radar  problem,  a 
procedure  similar  to  the  CFAR  approach  could  be  utilized  in  which  this  reference 
data  is  collected  from  range  cells  adjacent  to  the  test  cell.  Preliminary 
consideration  was  given  to  this  procedure  in  the  single  channel  case  [14,  16]  with 
considerable  detection  performance  improvement  noted. 

Another  area  of  investigation  would  involve  model-fitting  of  vector 
observation  processes  with  multichannel  time  series  for  specific  applications.  For 
the  radar  case,  single  channel  clutter  processes  have  been  considered  in  terms  of 
autoregressive  processes  [20,  21]  The  extension  of  these  analyses  to  multichannel 
processes  holds  potential  for  further  research  investigations  which  may  provide 
performance  improvement.  In  addition,  the  multichannel  model-based  approach 
appears  to  offer  the  potential  to  utilize  data  from  distinct,  yet  correlated, 
processes  in  the  detection  problem;  ie.  if  a  vector  time  series  exists  which  models 
the  individual  processes,  the  likelihood  ratio  presented  here  should  be  applicable. 
The  above  extensions  will  be  investigated  as  time  permits. 
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APPENDIX  A 

In  this  section,  we  determine  the  conditions  under  which  £(n)  as  expressed 
in  eq  (3-1)  is  a  temporal  white  noise  MMSE  output  of  a  linear  prediction  error 
filter.  We  first  follow  an  argument  similar  to  the  abbreviated  discussion  in  [9]. 
The  linear  prediction  error  of  *(n)  as  defined  in  eq  (4-1)  is  expressed  as 

£(n)=*(n)-£(n|n-l)  (A-l) 

where  x(n|n-l)  represents  the  estimated  vector  of  &(n)  using  past  data  values  with 
the  initial  condition  &(1|0)  =  Q.  Using  a  linear  predictor  with  P  past  values,  we 
define 

P  H 

£(n|n-l)  =  -  X  AP(k)i(n-k)  (A-2) 

k=l 

H 

where  Ap(k),  k  =  1,2,...,P  are  JxJ  matrices  representing  the  coefficients  of  the 

linear  predictor.  Substituting  eq  (A-2)  in  (A-l) 

P  H 

£(n)  =  i(n)  +  X  Ap(k)i(n*k)  (A-3a) 

k=l 

P  H 

=  s  Ap(k)i(n-k)  (A-3b) 

k=0 

JJ 

where  Ap(0)  =  I.  Let  the  concatenated  column  vector  of  P+1  vectors  (each  of 
dimension  J)  be  defined  as 

*n,n-p  =  UT(n)  *T(n-l)...&T(n-P).  (A-4) 

J| 

Post  multiplying  eq  (A-3)  by  xn,n-P  and  taking  the  expected  value 

H  H  H 

E[£(n)£n,n-p]  =  E[Ap  4n,n-PAn,n-p] 

•AHrr  H  i 

=  ApE[4n>n_PAp>n.p] 
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=Ap[ft,x) 


(A-5) 


where  Ap  and  [  Rxxl  are  defined  in  eq  (4-3b)  and  (4-3c).  The  LHS  of  eq  (A-5) 

T 

can  also  be  written  using  the  Hermitian  of  Xn,n*P  s0  that 

E  [£  (n) * “n_p3  =  E  {£  (n)  [*%)  x%-l) ...  £%-P))} 

=  {E  fe(n)iH(n)]  E  [£(n)*H(n-l)] ...  E[£(nUH(n-P)]}. 

(A-6) 

H 

We  now  determine  the  coefficients  of  Ap  subject  to  the  condition  that  £(n)  is  a 


MMSE  vector.  Under  this  condition, 

E  fe  (n)  AH(n-k)]  =  [0]  k>0. 


(A-7) 


Eq  (A-7)  is  the  orthogonality  condition  which  states  that  the  error  vector  is 
orthogonal  to  past  observation  values.  Using  this  condition,  eq  (A-6)  becomes 


Efe(n)i”n_p]=<E  fe(n)sH(n)]  [0]  [0] ...  [0]}. 


(A-8) 


From  eq  (A-3a) 

P*  H 

i(n)  =  £(n)  -  X  AP(k)i(n-k) 
k=l 

so  that 

P 

*H(n)  =  £H(n)  -  X  XH(n-k)AP(k). 
k=l 

Using  eq  (A- 10)  in  the  RHS  of  eq  (A-8) 

H  r  p 

E[£(nU" n-pJ  =  E[£(n)£H(n)]-  X  E[£(n)*H(n-k)]AP(k) 

.  k=l 


(A-9) 


(A- 10) 


[0]...[0] 
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={[£f]e  [0]  [0]...[0]} 


(A-ll) 


where  we  have  again  used  eq  (A-7)  to  yield  eq  (A-ll)  and  (Zf]e  is  the  forward 
prediction  error  covariance  matrix.  Combining  eq  (A-5)  and  (A-ll),  we  have 

Ap(ftxi]  =  {[lf]e[0][0]...[0]}.  (A- 12) 

Eq  (A-12)  is  the  multichannel  AR  Yule  Walker  normal  equation  in 
augmented  form.  It  provides  a  set  of  JP  linear  equations  to  solve  for  the  values  of 
the  matrix  coefficients  which  minimize  the  mean  square  error  vector.  Although 
eq  (A-12)  has  often  been  presented  in  the  literature,  the  reversed  order  form  of 
the  correlation  matrix  has  not  often  been  sighted  [1].  We  will  utilize  this  feature 
in  Appendix  B. 

We  now  show  that  the  vector  process  £(n)  is  uncorrelated  in  time.  At  an 
arbitrary  time  (n-l)  where  l>0,  eq  (A- 10)  becomes 

P 

*H(n-l)  =  £H(n-l)  -  X  *H(n-k-l)Ap(k).  I  >  0.  (A-13) 

k=l 

Using  eq  (A-13)  in  eq  (A-7) 

p 

E[£(n)aH(n-l)]  =  E[£(n)£H(n-l)]-  £  E[£(i.)AH(n-k-l)]AP(k)  =  [0].  (A-14) 

k=l 

From  eq  (A-7),  we  have 

E  [£  (n)  i  H  (n-k-l)]  =  [0]  k  >  0 .  ( A,  j  ^ 

And  so,  eq  (A-14)  yields 

E  [£  (n)£H(n-0  =  [0]  l>0.  (A-16) 

Thus,  the  sequence  of  outputs  from  the  MMSE  prediction  error  filter  are 
orthogonal.  Since  £(n)  is  a  zero  mean  Gaussian  process,  its  sequence  of  values 
{£(n)}  are  mutually  independent  so  that  {£(n)}  is  a  white,  Gaussian  noise 
sequence. 
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APPENDIX  B 


In  this  Appendix,  we  show  that  the  matrix  coefficients  of  multichannel 
linear  prediction  for  a  multichannel  random  process  and  the  prediction  error 
covariance  matrices  are  related  to  the  covariance  matrix  through  a  block 
triangular  decomposition.  The  procedure  is  a  straight  forward  generalization  of 
[10].  In  this  discussion,  we  are  considering  the  first  stage  of  processing  as  noted 
in  sections  IV  through  VI.  This  stage  results  in  the  £(n)  vector  output.  Therefore, 

we  only  consider  the  the  A(k)  matrix  coefficients  here.  A  treatment  which 
includes  the  second  stage  of  processing  using  the  matrix  Ly  to  obtain  y(n)  is 

developed  in  [36]. 

Recognizing  that  defined  in  eq  (2-5)  is  Hermitian  with  non-singular 
upper  left  principal  minors,  we  can  obtain 


db  _  i  n  i  H 

^XX-  £ 


(B-l) 


where  is  lower  block  triangular  with  the  identity  matrix  I  forming  the  block 
diagonal  matrices  and  Dg  is  a  real  block  diagonal  matrix.  Solving  eq  (B-l)  for 
D£ 


-1  B  ,H 
Dfi-L^U1  ). 

If  we  consider 

£tfT  *  i,n 
we  can  easily  show  that 
DP  =  E  [e  eh  ] 

£  “  1,N  “  1,N  J 

where 


(B-2) 


(B-3) 


(B-4) 


56 


£1^=  [£  0)£  (2)  —£  (N)  ]T 


(B-5) 


and  £(k)  is  a  J  x  1  channel  vector.  Since  Lg  is  lower  block  triangular  with  unit 
diagonal  elements,  Lg-1  has  the  same  form  so  that  eq  (B-3)  is  a  causal  and 
causally  invertible  transformation  of  the  data. 

We  now  consider  the  normal  equations  for  a  multichannel  predictor  of 
order  p  such  that 

A?[ftii]  =  (Ef]?IT  (B-6) 


where  [  R^lp  is  the  reversed  order  multichannel  covariance  matrix, 

A?  =  [IA?(1)  Ap(2)...Ap(P)J 


(B-7) 


IT  =  {I[0][0]...[0]}  (B-8) 

J.J 

where  I  is  a  JxJ  identity  matrix.  The  vector  of  matrices  Ap  is  the  vector  of 

multichannel  pth  order  linear  prediction  coefficients  and  [Lfjp  is  the  Hermitian, 

pth  order,  multichannel  prediction  error  covariance  matrix.  Post  multiplying  eq 
(B-6)  by  AP,  and  recognizing  that  [£f]p  is  Hermitian,  we  obtain 

[IA?(1)  A?(2)...A?(P)]fRJJ'  I  1 


Ap(l) 

Ap(2) 


AP(P) 


=  [Zf]P. 


(B-9) 


Using  the  relation 


G'1  G  =  I 

B  B 


(B-10) 


where  Gb  is  the  block  reflection  matrix  (i.e.,  the  square  J(p+1)  x  J(p+1)  matrix 
with  JxJ  identity  matrices  along  the  block  cross  diagonal),  one  can  easily  show 
that 


57 


Ap(P) 


[A?(P)...A?(2)  Ap(J)I][RMJ 


Ap(2) 

AP(1) 

-  I  - 


=  ffflp. 


(B-l  1) 


At  this  point,  it  is  noted  that  eq  (B-l  1 )  results  because  of  the  reversed 
order  of  the  correlation  matrix  [1]. 

We  now  write  eq  (B-l  1)  for  p  =  0,  1,  ...,  P  =  N-l  in  the  combined  form  as 


I 

R(0)  R(-l)  R(-2) ...  R(-P) 

a"(1)  I  t°l 

R(l)  R(0)  R(l) ...  R(-P+l) 

a" (2)  A?(l)  I 

R(2)  R(l)  R(0)  ...R(-P+2) 

a"(3)  A"  (2)  A?(l)  I  ... 

• 

• 

: 

• 

Ap(P)  A?(P-1)  a£(P-2)...  j$(1)J 

R(P)  R(P-l)  R(P-2)...R(0) 

IA,(1)  A^2)  Afi)  ...A^P) 

^1 

<-A 

o 

© 

I  Aj(l)  A^2)  ...Ap(P-l) 

fff], 

I  AJ1)  ...AP(P-2) 

Bf  12 

(0)  i  ..  : 

'  I 

_  101  BfIp 

(B-12) 
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Eq  (B-12)  is  of  the  same  form  as  eq  (B-2).  Since  the  causal  decomposition 
in  eq  (B  -2)  is  unique* ,  then  Lg'l  can  be  identified  with  the  lower  triangular 
matrix  in  eq  (B-12).  Thus  the  block  rows  of  are  the  multichannel 
coefficients  for  linear  predictive  filters  of  orders  zero  through  P  =  N-l  and  the 
block  diagonal  matrix  elements  of  are  the  prediction  error  variances 
associated  with  the  multichannel  filter  orders. 


*  The  uniqueness  of  this  decomposition  is  based  upon  the  specified  block  form  of  R®Xx  defined  in 
eq  (2-5).  However,  other  block  forms  of  could  have  been  made  which  still  retain  the 
Hermitian  property. 
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APPENDIX  C 


In  this  Appendix,  a  recursive  expression  for  the  sample  variance  estimate 
of  the  complex  quantity  Yj(n|Hj)  is  derived.  The  sample  mean  of  the  jth  channel 

sequence  Yj(n|Hi)  is  given  at  the  Nth  time  sample  as 

N 

7j(N|Hi)=  -L-  £  YjfllHj)  1=0,1.  (C-l) 

N  I  =1 

By  definition,  the  sample  variance  of  the  complex  quantity  ^(n|Hi)  is  expressed  as 

N  2 

$j(N|Hj)  =  j±-  V  lYj(klHi)  -Vj(NlHj)  |  1=0,1  (C-2a) 


N 

J_  f 
N-l  kTi 


N 

y 

N-lfl 


[  Yj(klHi)  -YjfNlHOl  lYfW^(NIHi)]  i=0,l 

(C-2b) 

(iTjCklHj)  l^jCklHj)  7*(N|Hi)-Yj(N|Hi)Yf(k|Hi) 

+|Yj(N|Hi)|2}  i-0,1. 

(C-2c) 


Using  eq  (C-l)  in  (C-2c) 
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^2(N|H')  =  Fn  |  |  Tj«H) 


-  Tr|  1'j(,|Hi)(FrT)k|  TjWV 

+  i^i  £  "n2  £  VjOIH,)!  y*(m|Hj)  i=0,l  (C-3) 

k=l  1=1  m=l 

-AS^i-sfe,  £  £  Yi(k,Hi)if?"H D 

1  V1  ^  1  N,  N,  N, 

•  N(nT)  2.  L  Y j 0 1 Hj )  Y *|(k|H j)  +  X  L  ZiVj(l|H,)Y*(m|H1) 

i=0,l . 

(C-4) 

Interchanging  k  and  I  in  the  second  summation  as  well  as  k  and  m  in  the  fourth 
does  not  change  the  expression.  And  so, 

£  |  Yj  (k|H.)  £  YjC  |Hj)Yj(k|Hj) 

+  £,  N^N-l)  XX  YjC |Hj)  Y*(k|Hj)  i-0,1. 

(C-5) 


But  the  last  summation  in  eq  (C-5)  is  just 
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so  that 


N  N 


jrfc)§gWY;«Hi) 


N  N 


=  _L_)^|Y.(k|Hl)|  .  _I__  £  £  YjfllHj)  Yj(k|Hj) 


N(N-l)  £ 


1  1=1 


i=0,l.  (C-6) 


To  obtain  a  recursive  form  for  aj2(N|Hi)  we  separate  out  the  terns  involving 
Yj(N|Hi)  so  that 


"n<nT)  ^k|H‘> 


N-l  N-l 


i^irJtNIHp  X  Iy/uh,)  t;«h,> 


i=0,l  (C-7a) 
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=(FM)lYi(NiHi)|2+  FTT  S  |Tj«HJ  f-  |7j(N|Hj)| 


fen  yj(N|Hi)  -^y^h^oih,) 


N(N-l) 


1  SW 

N(N-l)  21  ^j(MHi)  YjCklHj)  i=0,l 


(C-7b) 


so  that 


A2(N|H)=  JL  |7j(N|H,)|  +  J_  |Yj(k|Hj)  |2-  J-T-.tNIHi)  Y/N-1|H,) 

ic—  1 

-TT^NIH.)  7j(N-l|Hj) -JjjfejjX  2 Yj« IHj) Y*(k|Hp  1-0,1 

I  1  k—  1 


(C-7c) 


From  eq  (C-6) 


’j!(N-"^,=(NL2,|!|7j«H^|2-  S  |lYj(MHi)Y;(k|Hi) 


i=0,l. 

(C-8) 


Multiplying  both  sides  by  N-2/N-1, 
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N-l  (N-l)^lYj(klH*) !  "(N-l)^Yj(l,Hi)(N-l)^  Yj(klHp 


i=0,l 

(C-9) 


so  that 


j~7  Z  iTjdcIHj)  f-  Oj2(N-l|H)  -^(N-llHj)  ?*(N-1|H,)  ,=0,1 


^■c/(N-'|Hi)+|7j(N-l|Hi)|2  U0,1, 


(C-10a) 


(C-lOb) 


Using  eq  (C-lOb)  in  eq  (C-7c) 


^(N!H)=Jr  |t'j(N|H,)|  aj2(N..m)+|yj(N-l|H,)| 

Y^NIHj)  Y^N-HHj)  TjtNIHj)  T^N-llHj) 

N  N 


(N-l) 

(N) 


[(N-l)  ^jYJ(l|Hi^  [(N-l)  S  YjOcIFii)] 


i=0,l 


(C-ll) 
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»r" 


=  F)  'Yj(NIHj)l2  +  ^f^(N-IIHj)  +  ^(N-llHj)!2 
-j5j  Yj(NIHj)  Y*(N-llHj)  +  Y* (NIHj)  Yj(N-llHj)  -  ^  ^(N-lIH;)!2 

so  that 

Oj2(NIHi)  =  ^o'VlIHj)  +  ^  lyj(NIHi)l2 

-  g  [Y/NIHi)?*(N-llHi)  +  T*(NIHi)Yj(N-llHj)| 

+  FjlYj(N-lllHi)l2  (C-I2a) 

=  ^^(N-llHj)  +  tYj(NIHi)  -  7j(N-llHil2  (C-12b) 

with  initial  conditions 

d?(llHi)  =  0  i  =  0,1.  (C-13) 
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